d03edf

d03edf © Numerical Algorithms Group, 2002.

Purpose

D03EDF Elliptic PDE, solution of finite difference equations by a multigrid technique

Synopsis

[a,rhs,us,u,numit,ub,ifail] = d03edf(ngx,ngy,a,rhs<,ub,maxit,acc,iout,ifail>)

Description

 
 D03EDF solves, by multigrid iteration, the seven-point scheme
 
  6            7                
 A   u       +A   u             
  i,j i-1,j+1  i,j i,j+1     
 
   3           4        5       
 +A   u      +A   u   +A   u    
   i,j i-1,j   i,j ij   i,j i+1,j
 
               1          2       
             +A   u     +A   u        =f  ,
               i,j i,j-1  i,j i+1,j-1  ij
 
                                      i=1,2,...,n ; j=1,2,...,n ,
                                                 x             y
 
 which arises from the discretization of an elliptic partial 
 differential equation of the form
 
 
 (alpha)(x,y)U  +(beta)(x,y)U  +(gamma)(x,y)U  +(delta)(x,y)U 
              xx             xy              yy              x
 
 +(epsilon)(x,y)U +(phi)(x,y)U=(psi)(x,y)
                 y
 
 and its boundary conditions, defined on a rectangular region. 
 This we write in matrix form as
 
                               Au=f
 
 Systems of linear equations, matching the seven-point stencil 
 defined above, are solved by a multigrid iteration. An initial 
 estimate of the solution must be provided by the user. A zero 
 guess may be supplied if no better approximation is available.
 
 A 'smoother' based on incomplete Crout decomposition is used to 
 eliminate the high frequency components of the error. A 
 restriction operator is then used to map the system on to a 
 sequence of coarser grids. The errors are then smoothed and 
 prolongated (mapped onto successively finer grids). When the 
 finest cycle is reached, the approximation to the solution is 
 corrected. The cycle is repeated for MAXIT iterations or until 
 the required accuracy, ACC, is reached.
 
 D03EDF will automatically determine the number l of possible 
 coarse grids, 'levels' of the multigrid scheme, for a particular 
 problem. In other words, D03EDF determines the maximum integer l 
 so that n  and n  can be expressed in the form
          x      y                         
 
               l-1          l-1                   
          n =m2   +1,  n =n2   +1,  with  m>=2 and n>=2.
           x            y                         
 
 It should be noted that the rate of convergence improves 
 significantly with the number of levels used, so that n  and n  
                                                        x      y
 should be carefully chosen so that n -1 and n -1 have factors of 
                                     x        y
           1
 the form 2 , with l as large as possible. For good convergence 
 the integer l should be at least 2.
 
 D03EDF has been found to be robust in application, but being an 
 iterative method the problem of divergence can arise. For a 
 strictly diagonally dominant matrix A
 
                          4    --    k
                        |A  |> >   |A  |
                          ij   --    ij
                               k/=4

 no such problem is foreseen. The diagonal dominance of A is not a
 necessary condition, but should this condition be strongly 
 violated then divergence may occur. The quickest test is to try 
 the routine.
 

Parameters

d03edf

Required Input Arguments:

ngx                                   integer
ngy                                   integer
a (:,7)                               real
rhs (:)                               real

Optional Input Arguments:                       <Default>

ub (ngx*ngy)                          real     zeros(ngx*ngy,1)
maxit                                 integer  15
acc                                   real     sqrt(eps)
iout                                  integer  0
ifail                                 integer  -1

Output Arguments:

a (:,7)                               real
rhs (:)                               real
us (:)                                real
u (:)                                 real
numit                                 integer
ub (ngx*ngy)                          real
ifail                                 integer